Code to load tools and data

  # TOOLS 
    require(here)
    source(here::here('R/tools.R'))
    
    require(arm) 
    require(effects)
    require(ggpubr)
    require(ggsci) 
    require(grid)
    require(gridExtra)
    require(magrittr)
    require(multcomp)
    require(PerformanceAnalytics)
    require(rptR) 
    require(stringi)
    require(viridis)

    # constants
      round_ = 3 # number of decimal places to round model coefficients
      nsim = 5000 # number of simulations to extract estimates and 95%CrI
      ax_lines = "grey60" # defines color of the axis lines
      colors <- c("#999999", "#E69F00", "#56B4E9") #viridis(3)
    
    # functions
      getime = function (x) {ifelse(is.na(x), as.numeric(NA), as.numeric(difftime(x, trunc(x,"day"), units = "hours")))}
      
      getDay = function (x) {as.Date(trunc(x, "day"))}
     
      # for R output based on sim
        R_out = function(name = "define", model = m, nsim = 5000){
         bsim <- sim(model, n.sim=nsim)  
         l=data.frame(summary(model)$varcor)
         l = l[is.na(l$var2),]
         l$var1 = ifelse(is.na(l$var1),"",l$var1)
         l$pred = paste(l$grp,l$var1)

         q50={}
         q025={}
         q975={}
         pred={}
         
         # variance of random effects
         for (ran in names(bsim@ranef)) {
           #ran =names(bsim@ranef)[1]
           ran_type = l$var1[l$grp == ran]
           for(i in ran_type){
              # i = ran_type[2]
            q50=c(q50,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.5)))
            q025=c(q025,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.025)))
            q975=c(q975,quantile(apply(bsim@ranef[[ran]][,,i], 1, var), prob=c(0.975)))
            pred= c(pred,paste(ran, i))
            }
           }
         # residual variance
         q50=c(q50,quantile(bsim@sigma^2, prob=c(0.5)))
         q025=c(q025,quantile(bsim@sigma^2, prob=c(0.025)))
         q975=c(q975,quantile(bsim@sigma^2, prob=c(0.975)))
         pred= c(pred,'Residual')

         ci = c(round(100*q025/sum(q025))[1], round(100*q975/sum(q975))[1])
         ci = ci[order(ci)]
         
         ri=data.table(model = name, repeatability=paste0(round(100*q50/sum(q50)),'%')[1], CI = paste0(paste(ci[1], ci[2], sep ="-"), '%'))
         
         
         return(ri)
         }
  # DATA 
    # composite measures
      x = fread(here::here('Data/DAT_morpho.csv')) 
      setnames(x,old = 'pic', new = 'sperm_ID')
      x[, sample_ID:=as.character(sample_ID)]
      
      # for each bird 10 sperm measured from one of the two sampling occasions
        #bb = b[part == 'Tail']
        #bbx = data.table(table(bb$bird_ID,bb$month))
        #bbx[!N%in%c(0,10), unique(V1)]

        #  mneasurements from minip vs rest are the same - so no need to control for
        #ggplot(x[part == 'Tail'], aes(x = manip, y = Pixels)) + geom_boxplot()
      
      d1 = x[,.(bird_ID, sample_ID, sperm_ID, part, Pixels)]
      
      d2 = d1[, sum(Pixels), by = list(bird_ID, sample_ID, sperm_ID)]
      names(d2)[4] = 'Pixels'
      d2[,part := 'Total']

      d3 = d1[part%in%c('Midpiece', 'Tail'), sum(Pixels), by = list(bird_ID, sample_ID, sperm_ID)]
      names(d3)[4] = 'Pixels'
      d3[,part := 'Flagellum']

      d4 = d1[part%in%c('Acrosome', 'Nucleus'), sum(Pixels), by = list(bird_ID, sample_ID, sperm_ID)]
      names(d4)[4] = 'Pixels'
      d4[,part := 'Head']

      b = merge(d1,d2, all = TRUE)
      b = merge(b,d3, all = TRUE)
      b = merge(b,d4, all = TRUE)
      b[, Length_µm:=Pixels*0.078]

    # add metadata 
      b = merge(b, x[,.(bird_ID, sample_ID, sperm_ID, part, manip)], all.x = TRUE, by=c('bird_ID', 'sample_ID', 'sperm_ID', 'part'))
      b[manip%in%"", manip:=NA]
      # adjust multiple samples from the same time to the one recorded for motility
        b[sample_ID==51, sample_ID:=52]
        b[sample_ID==3, sample_ID:=4]
      s = data.table(read_excel(here::here('Data/sampling_2021_cleaned.xlsx'), sheet = 1))
      s = s[!is.na(recording)]
      m = data.table(read_excel(here::here('Data/ruff_males_Seewiesen.xlsx'), sheet = 1))#, range = "A1:G161"))
      m[, hatch_year:=as.numeric(substr(Hatch_date,1,4)) ]
      m[, age := 2021-hatch_year]

      s = merge(s,m[,.(Ind_ID, Morph, age)], by.x = 'DB_ID', by.y = 'Ind_ID', all.x = TRUE)
      # adjust IDs (where missed typed) for smooth merging
        b[bird_ID=='A027121649', bird_ID:='AO27121649']
        b[bird_ID=='A03999183', bird_ID:='AO3999183']
        b[bird_ID=='A0414173NL', bird_ID:='AO414-17-3NL']
        b[bird_ID=='A0414175', bird_ID:='AO414-17-5']
        b[bird_ID=='A059381830', bird_ID:='AO59381830']
        b[bird_ID=='A07422181', bird_ID:='AO7422181NL']
        b[bird_ID=='A079561654', bird_ID:='AO79561654']
        b[bird_ID=='A079561656', bird_ID:='AO79561656']
        b[bird_ID=='A079561718', bird_ID:='AO7956-17-18']
        b[bird_ID=='AIFA-016507', bird_ID:='AIFAO16507']
        b[bird_ID=='AIFA01511', bird_ID:='AIFAO15-11']
        b[bird_ID=='A079561718', bird_ID:='AO79561718']
        b[bird_ID=='G20005', bird_ID:='G200055']
        b[bird_ID=='cz005', bird_ID:='CZ005']
        b[bird_ID=='AO414-17-3NL', bird_ID:='AO414-17-3-NL']

      b = merge(b, s[,.(sample_ID, bird_ID, Morph, age, datetime, type, sperm, recording, rec_measured, month)], by.x = c('sample_ID', 'bird_ID'), by.y = c('sample_ID', 'bird_ID'), all.x = TRUE)
      
      # fwrite(ss[!duplicated(bird_ID), .(DB_ID, bird_ID)], file = 'Data/used_males_for_Clemens.csv')
      #b = (b[!duplicated(b)]) # shouldn't be necessary,
    
    # add motility on to individual morpho measurements (only one motility value per individual)
      d = data.table(read_excel(here::here('Data/motility.xlsx'), sheet = 1))
      d[, motileCount_ln := log(motileCount)]
      setnames(d, old=c('date','ID'), new = c('month','bird_ID'))

      # for bird 1339 motility measured from May sample, but sperm from June, so meta data changed to May sample 
         b[sample_ID==183, month:='May']
         b[sample_ID==183, datetime:=as.POSIXct('2021-05-04 12:33')]
         b[sample_ID==183, sample_ID:=95]
      
      b = merge(b, d, by = c('bird_ID', 'month'), all.x = TRUE)

      #b[is.na(Morph), unique(bird_ID)]
     
      b[Morph == 'F', Morph := 'Faeder']
      b[Morph == 'I', Morph := 'Independent']
      b[Morph == 'S', Morph := 'Satellite']
      b[is.na(issues), issues := 'zero']

      b[, part := factor(part, levels = c('Acrosome', 'Nucleus', 'Midpiece', 'Tail','Head','Flagellum', 'Total'))]
      b[part %in% c('Acrosome', 'Nucleus', 'Midpiece', 'Tail'), measure := 'original']
      b[part %in% c('Head','Flagellum', 'Total'), measure := 'composite']

      #nrow(d) # N 139
   
    # add metadata to motility dataset
      # check for multiple samples per individuals at a single sampling period & keep only one for merging
        ajm = data.table(table(s$bird_ID,s$month))
        #s[bird_ID%in%ajm[N>1, unique(V1)]]
        ss = s[rec_measured %in% 'yes']    
      d = merge(d, ss[,.(sample_ID, bird_ID, Morph, age, datetime, type, sperm, recording, month)], by = c('month','bird_ID'), all.x = TRUE)
     
      d[is.na(Morph), Morph := 'Zebra finch']
      d[Morph == 'F', Morph := 'Faeder']
      d[Morph == 'I', Morph := 'Independent']
      d[Morph == 'S', Morph := 'Satellite']
      d[is.na(issues), issues := 'zero']

    # prepare for correlations and repeatability
      bw = reshape(b[,.(bird_ID,Morph, age, datetime, month, sample_ID, sperm_ID, VAP,VSL,VCL, motileCount, part, Length_µm)], idvar = c('bird_ID','Morph','age','datetime', 'month', 'sample_ID', 'sperm_ID','VAP','VSL','VCL', 'motileCount'), timevar = 'part', direction = "wide")  
      setnames(bw,old = c('Length_µm.Acrosome', 'Length_µm.Nucleus','Length_µm.Head','Length_µm.Midpiece','Length_µm.Tail','Length_µm.Flagellum', 'Length_µm.Total'), new = c('Acrosome', 'Nucleus', 'Head','Midpiece', 'Tail','Flagellum','Total'))
      # add relative measures
        bw[, Midpiece_rel := Midpiece/Total]
        bw[, Flagellum_rel := Flagellum/Total]

     dw = reshape(d[bird_ID%in%d[duplicated(bird_ID), bird_ID],.(bird_ID,species, Morph, age, month, VAP,VSL,VCL, motileCount, motileCount_ln)], idvar = c('bird_ID','species','Morph','age'), timevar = 'month', direction = "wide")  
   
    # mean/male dataset
      a = b[, list(mean(Length_µm), mean(VAP),mean(VSL), mean(VCL), mean(motileCount)), by = list(month, bird_ID, Morph, age, part)]
       setnames(a, old = c('V1', 'V2','V3','V4','V5'), new = c('Length_avg', 'VAP', 'VSL', 'VCL', 'motileCount'))
      a1 = data.table(bird_ID = unique(a$bird_ID), blind = c(rep(c('Independent','Satellite', 'Faeder'), floor(length(unique(a$bird_ID))/3)), 'Independent','Satellite'))
      a = merge(a,a1, all.x = TRUE)

      aw = reshape(a, idvar = c('month','bird_ID','Morph', 'blind','age','VAP','VSL','VCL', 'motileCount'), timevar = "part", direction = "wide")
      names(aw) = c('month','bird_ID','Morph', 'VAP','VSL','VCL', 'motileCount', 'blind','age',as.character(unique(a$part)))
      # add relative measures
      aw[, Midpiece_rel := Midpiece/Total]
      aw[, Flagellum_rel := Flagellum/Total]

    # morph as factor
      b[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))] 
      bw[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))] 
      a[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))] 
      aw[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder"))] 
      
      d[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder", "Zebra finch"))] 
      dw[, Morph := factor(Morph, levels=c("Independent", "Satellite", "Faeder", "Zebra finch"))] 
      #ggplot(dw, aes(x = VAP.May, y = VAP.June))+geom_point() + stat_smooth(method = 'lm') + geom_abline(intercept = 0, slope = 1, lty = 3, col = 'red')
    # dataset for relative measurements
      bt = b[part == 'Total',.(Morph, bird_ID, sample_ID, sperm_ID, part, Length_µm, VAP, VSL, VCL, motileCount, motileCount_ln)]
      
      b1 = b[part%in%c('Midpiece'),.(Morph, bird_ID, sample_ID, sperm_ID, part, Length_µm, VAP, VSL, VCL, motileCount, motileCount_ln)]
      b1[, Length_rel := Length_µm/bt$Length_µm]
     
      b2 = b[part%in%c('Flagellum'),.(Morph, bird_ID, sample_ID, sperm_ID, part, Length_µm, VAP, VSL, VCL, motileCount, motileCount_ln)]
      b2[, Length_rel := Length_µm/bt$Length_µm]
      
      br = rbind(b1,b2)
      br$Length_µm = NULL

      at = a[part == 'Total',.(Morph, bird_ID,  part, Length_avg, VAP, VSL, VCL, motileCount)]
      
      a1 = a[part%in%c('Midpiece'),.(Morph, bird_ID, part, Length_avg, VAP, VSL, VCL, motileCount)]
      a1[, Length_rel := Length_avg/at$Length_avg]
     
      a2 = a[part%in%c('Flagellum'),.(Morph, bird_ID, part, Length_avg, VAP, VSL, VCL, motileCount)]
      a2[, Length_rel := Length_avg/at$Length_avg]
      
      ar = rbind(a1,a2)
      ar$Length_avg = NULL
    
    # dataset for correlations
      h = b[part == 'Head']
      setnames(h, old = "Length_µm", new="Head_µm")
      h[, Acrosome_µm := b[part == 'Acrosome',.(Length_µm)]]
      h[, Nucleus_µm := b[part == 'Nucleus',.(Length_µm)]]
      h[, Midpiece_µm := b[part == 'Midpiece',.(Length_µm)]]
      h[, Tail_µm := b[part == 'Tail',.(Length_µm)]]
      h[, Flagellum_µm := b[part == 'Flagellum',.(Length_µm)]]
      h[, Total_µm := b[part == 'Total',.(Length_µm)]]
      h[, Midpiece_rel := br[part == 'Midpiece',.(Length_rel)]]
      h[, Flagellum_rel := br[part == 'Flagellum',.(Length_rel)]]

      ha = a[part == 'Head']
      setnames(ha, old = "Length_avg", new="Head_µm")
      ha[, Acrosome_µm := a[part == 'Acrosome',.(Length_avg)]]
      ha[, Nucleus_µm := a[part == 'Nucleus',.(Length_avg)]]
      ha[, Midpiece_µm := a[part == 'Midpiece',.(Length_avg)]]
      ha[, Tail_µm := a[part == 'Tail',.(Length_avg)]]
      ha[, Flagellum_µm := a[part == 'Flagellum',.(Length_avg)]]
      ha[, Total_µm := a[part == 'Total',.(Length_avg)]]
      ha[, Midpiece_rel := ar[part == 'Midpiece',.(Length_rel)]]
      ha[, Flagellum_rel := ar[part == 'Flagellum',.(Length_rel)]]
    
    # CV dataset
        cv_ =  b[, cv(Length_µm), by = list(bird_ID, part, Morph)]
        cv_[ , Morph123 := as.numeric(Morph)]
        names(cv_) [4]='CV'

Basic info

  • for each male 10 sperm measured from one of the two sampling occasions
  • motility for most individuals and for many in May and June
  • for all morpho data we have a motility measurement that corresponds with the sampling date, except for one male, where motility measured in May, but morpho from June sample
  • initial exploration of the motility is available from here.
  • to decide:
    • whether to control for pedigree
    • for motility ~ morph analyses I, use
      1. June motility values and for 4 males without June, May values. Is this ok or better to use male averages?
      2. bird_ID as random intercept - 50 males measured once 42 twice (June & May)
    • motility models are controlled for number of tracked sperm, do we need to control also for issues (like few sperm, shit, running sample) or age?
    • in motility ~ morphology - the motility measurement (with one exception) comes from the same sample as the morphology measurement. Is this ok?

Exploration

  b[, order_ := mean(Length_µm), by = bird_ID]
  b_ = b[part =='Total']
  b_[,bird_ID := reorder(bird_ID, Length_µm, mean)]
  b[, bird_ID := factor(bird_ID, levels = levels(b_$bird_ID))]

  g = ggplot(b, aes(x = reorder(as.factor(bird_ID),Length_µm,mean), y = Length_µm)) + facet_wrap(~part, nrow = 7, scales = 'free') + #x = reorder(as.factor(bird_ID),Length_µm,mean)
    geom_boxplot(aes(col = Morph)) + 
    labs(subtitle = 'Distribution within & across males (ordered by total length)')+
    #geom_dotplot(binaxis = 'y', stackdir = 'center',
     #            position = position_dodge(), col = 'red', fill ='red')+
    scale_colour_manual(values = colors) + 
    #scale_color_viridis(discrete=TRUE)+
    xlab('Male ID') +
    theme_bw() +
    theme(axis.text.x = element_blank())
    #$, legend.position = "none")
  g

  ggsave(here::here('Output/morpho_within_male_boxplots_ordered.png'), g, width = 20, height = 15, units = 'cm')
  #fig.width=5, fig.height = 5
  chart.Correlation(bw[, c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
  mtext("Single sperm", side=3, line=3)
  #dev.copy(png,'Output/VD_corr_single.png')
  #dev.off()

  # cor_parts_avg, fig.width=5, fig.height=5  
  chart.Correlation(aw[, c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total','Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
  mtext("Male averages", side=3, line=3)
  #dev.copy(png,'Output/VD_corr_avg.png')
  #dev.off()

!!! Nucleus strongly predict head and tail total sperm length. !!!


  #fig.width=5, fig.height = 5  
  chart.Correlation(bw[Morph == 'Independent', c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
  mtext("Independent", side=3, line=3)
  #dev.copy(png,'Output/morpho_corr_single_I.png')
  #dev.off()  
  #cor_parts_S, fig.width=5, fig.height = 5    
  chart.Correlation(bw[Morph == 'Satellite', c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
  mtext("Satellite", side=3, line=3)
  #dev.copy(png,'Output/morpho_corr_single_S.png')
  #dev.off()
  #cor_parts_F, fig.width=5, fig.height = 5    
  chart.Correlation(bw[Morph == 'Faeder', c('Acrosome', 'Nucleus','Head', 'Midpiece', 'Tail', 'Flagellum','Total', 'Midpiece_rel', 'Flagellum_rel')], histogram=TRUE, pch=19)
  mtext("Feader", side=3, line=3)
  #dev.copy(png,'Output/morpho_corr_single_F.png')
  #dev.off()

!!! Correlations look similar across morphs, except for midpiece in faeders where it more strongly correlates with tail and hence also total sperm length !!!


Repeatability

  # male 
    lfsim = list()
    lfrpt = list()
    for(i in c('Acrosome','Flagellum','Head','Midpiece','Nucleus','Tail','Total')){
      part_ = i
      # part_ = "Acrosome"
      dd = b[part == part_]
      m = lmer(Length_µm ~ 1+(1|bird_ID), dd)
      Rf = R_out(part_)
      lfsim[[i]] = Rf[, method_CI:='arm package']

      R = rpt(Length_µm ~ (1 | bird_ID), grname = "bird_ID", data = dd, datatype = "Gaussian")#, nboot = 0, npermut = 0)
      RR = data.table(merge(data.frame(name =part_), paste0(round(R$R*100),'%'))) %>% setnames(new = c('part', 'Repeatability'))
      RR[, CI := paste0(paste(round(R$CI_emp*100)[1], round(R$CI_emp*100)[2], sep = "-"), '%')] 
      lfrpt[[i]] =  RR[, method_CI := 'rpt package']
      print(i)
    }
    
    x = do.call(rbind,lfsim)
    names(x)[1] = "part"
    y = do.call(rbind,lfrpt)

    x[, pred:= as.numeric(substr(repeatability,1,2))]
    x[, lwr:= as.numeric(substr(CI,1,2))]
    x[, upr:= as.numeric(substr(CI,4,5))]

    y[, pred:= as.numeric(substr(Repeatability,1,2))]
    y[nchar(CI) == 5, CI := paste0(0,CI) ]
    y[, lwr:= as.numeric(substr(CI,1,2))]
    y[, upr:= as.numeric(substr(CI,4,5))]
    names(y)[2] = tolower( names(y)[2])
    xy = rbind(x,y)
    xy[, part := factor(part, levels=c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total"))] 
  # morph
    lfsim = list()
    lfrpt = list()
    for(i in c('Acrosome','Flagellum','Head','Midpiece','Nucleus','Tail','Total')){
      part_ = i
      # part_ = "Acrosome"
      dd = b[part == part_]
      m = lmer(Length_µm ~ 1+(1|Morph), dd)
      Rf = R_out(part_)
      lfsim[[i]] = Rf[, method_CI:='arm package']

      R = rpt(Length_µm ~ (1 | Morph), grname = "Morph", data = dd, datatype = "Gaussian")#, nboot = 0, npermut = 0)
      RR = data.table(merge(data.frame(name =part_), paste0(round(R$R*100),'%'))) %>% setnames(new = c('part', 'Repeatability'))
      RR[, CI := paste0(paste(round(R$CI_emp*100)[1], round(R$CI_emp*100)[2], sep = "-"), '%')] 
      lfrpt[[i]] =  RR[, method_CI := 'rpt package']
      print(i)
    }
    
    x = do.call(rbind,lfsim)
    names(x)[1] = "part"
    y = do.call(rbind,lfrpt)

    x[, pred:= as.numeric(sub("\\%.*", "", repeatability))]
    x[, lwr:= as.numeric(sub("\\-.*", "", CI))]
    x[, upr:= as.numeric(sub("\\%.*", "",sub(".*-", "", CI)))]
   
    y[, pred:= as.numeric(sub("\\%.*", "", Repeatability))]
    y[, lwr:= as.numeric(sub("\\-.*", "", CI))]
    y[, upr:= as.numeric(sub("\\%.*", "",sub(".*-", "", CI)))]
    #y[nchar(CI) == 5, CI := paste0(0,CI) ]
    #y[, lwr:= as.numeric(substr(CI,1,2))]
    #y[, upr:= as.numeric(substr(CI,4,5))]
    names(y)[2] = tolower( names(y)[2])
    xy2 = rbind(x,y)
    xy2[, part := factor(part, levels=c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total"))] 
  g1 = 
      ggplot(xy, aes(x = part, y = pred, col = method_CI)) +
        geom_errorbar(aes(ymin = lwr, ymax = upr, col = method_CI), width = 0.1, position = position_dodge(width = 0.25) ) +
        #ggtitle ("Sim based")+
        geom_point(position = position_dodge(width = 0.25)) +
        scale_color_viridis(discrete=TRUE, begin=0, end = 0.5)  +
        scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) + 
        labs(x = NULL, y = "Repeatability [%]", subtitle = "within male")+
        ylim(c(0,100))+
        coord_flip()+
        theme_bw() +
        theme(plot.title = element_text(size=9),
          legend.position = "none")

  #ggsave(here::here('Output/morpho_Repeatability_within-male.png'),g1, width = 10, height =7, units = 'cm')

  g2 = 
      ggplot(xy2, aes(x = part, y = pred, col = method_CI)) +
        geom_errorbar(aes(ymin = lwr, ymax = upr, col = method_CI), width = 0.1, position = position_dodge(width = 0.25) ) +
        #ggtitle ("Sim based")+
        geom_point(position = position_dodge(width = 0.25)) +
        scale_color_viridis(discrete=TRUE, begin=0, end = 0.5, guide = guide_legend(reverse = TRUE)) +
        #scale_fill_viridis(discrete=TRUE, begin=0, end = 0.5) + 
        labs(x = NULL, y = "Repeatability [%]", subtitle = "within morph")+
        ylim(c(0,100))+
        coord_flip()+
        theme_bw() +
        theme(plot.title = element_text(size=9),
          axis.title.y = element_blank(), 
        axis.text.y = element_blank())
    
    #ggsave(here::here('Output/morpho_Repeatability_within-morph.png'),g, width = 10, height =7, units = 'cm')

  grid.draw(cbind(ggplotGrob(g1), ggplotGrob(g2), size = "last"))

  #grid.arrange(g1,g2)
  ggsave(here::here('Output/morpho_Repeatability.png'),cbind(ggplotGrob(g1),ggplotGrob(g2) , size = "last"), width = 6*2, height =3*1.5, units = 'cm')  

Differences - “raw”

Red indicates mean

    g1 =  # dummy to extract variables for median calculation
      ggplot(b, aes(x = Morph, y = Length_µm)) +
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 0.5)+
      geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
      stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
      scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
      scale_color_viridis(discrete=TRUE, alpha = 0.8)+
      facet_wrap(~part, scales = 'free_y', nrow = 2)+
      ggtitle('Single measurements') +
      ylab('Length [µm]') +
      theme_bw() +
      guides(x =  guide_axis(angle = -45)) +
      theme(legend.position = "none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) #panel.spacing.y = unit(0, "mm")) #, 
 
    g2 =  
      ggplot(a, aes(x = Morph, y = Length_avg)) +
      geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 1)+
      #              position = position_dodge(), col = 'grey', aes(fill =Morph), dotsize = 1)+
      scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
      scale_color_viridis(discrete=TRUE, alpha = 0.8)+
      stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
      facet_wrap(~part, scales = 'free_y', nrow = 2)+
      ggtitle('Male means') +
      ylab('Length [µm]') +
      theme_bw() +
      guides(x =  guide_axis(angle = -45)) +
      theme(legend.position = "none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) #panel.spacing.y = unit(0, "mm")) #axis.title.x = element_blank(), axis.text.x = element_blank(), 
    
    cv_ =  b[, cv(Length_µm), by = list(bird_ID, part, Morph)]
    cv_[ , Morph123 := as.numeric(Morph)]
    names(cv_) [4]='CV' 
    
    g3 = 
     ggplot(cv_, aes(x =Morph , y = CV)) + 
     geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
     geom_dotplot(binaxis = 'y', stackdir = 'center',
                     position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 1)+
     stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
     scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
     scale_color_viridis(discrete=TRUE, alpha = 0.8)+
     facet_wrap(~part, nrow = 2, scales = "free_y") +
     guides(x =  guide_axis(angle = -45)) +
     ylab('Coefficient of variation') +
     ggtitle('Male values') +
     theme_bw() +
     theme(legend.position = "none",
          plot.title = element_text(size=9),
          axis.title.x = element_blank()
          )
    
    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    
    grid.draw(rbind(gg1, gg2, gg3, size = "last"))

      #grid.arrange(g1,g2)
 
    ggsave(here::here('Output/morpho_boxplots_.png'), rbind(gg1,gg2,gg3, size = "last"), width = 7*3, height =10*1.5, units = 'cm')  

    #ggsave(here::here('Output/morpho_per_male.png'),g, width = 10, height =10, units = 'cm')
   g=
   ggplot(cv_, aes(x =Morph , y = CV, col = part, fill = part)) + 
   #geom_dotplot(binaxis = 'y', stackdir = 'center', position = position_dodge(), col = 'grey',  dotsize = 0.5)+
   geom_boxplot(fill = NA) + 
   theme_bw() +
   theme(axis.title.x = element_blank())
   g

   ggsave(here::here('Output/CV_per_male_alternative.png'),g, width = 10, height =7, units = 'cm')
    g1 =  # dummy to extract variables for median calculation
      ggplot(br, aes(x = Morph, y = Length_rel)) +
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 0.5)+
      geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
      stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
      scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
      scale_color_viridis(discrete=TRUE, alpha = 0.8)+
      facet_wrap(~part, scales = 'free_y', nrow = 2)+
      ggtitle('Single measurements') +
      ylab('Relative to total sperm length') +
      theme_bw() +
      guides(x =  guide_axis(angle = -45)) +
      theme(legend.position = "none",
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) #panel.spacing.y = unit(0, "mm")) #, 
    g2 =  # dummy to extract variables for median calculation
      ggplot(ar, aes(x = Morph, y = Length_rel)) +
      geom_dotplot(binaxis = 'y', stackdir = 'center',
                   position = position_dodge(), aes(col = Morph, fill =Morph), dotsize = 0.5)+
      geom_boxplot(col = 'grey40', fill = NA, alpha = 0.2) + 
      stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
      scale_fill_viridis(discrete=TRUE, alpha = 0.4)+
      scale_color_viridis(discrete=TRUE, alpha = 0.8)+
      facet_wrap(~part, scales = 'free_y', nrow = 2)+
      ggtitle('Male means') +
      ylab('Relative to total sperm length') +
      theme_bw() +
      guides(x =  guide_axis(angle = -45)) +
      theme(legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        plot.title = element_text(size=9)
        ) 
    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    grid.draw(cbind(gg1, gg2, size = "last"))

      #grid.arrange(g1,g2)
      
    ggsave(here::here('Output/morpho_boxplots_rel.png'),cbind(gg1,gg2, size = "last"), width = 7*2, height =10*1.5, units = 'cm')    

Morpho-space

  # single sperm  
    bs = b[, quantile(Length_µm, prob = c(0.5)), by = list(part,Morph)]
    names(bs)[3] = 'median'
    bs$lwr = b[, quantile(Length_µm, prob = c(0.025)), by = list(part,Morph)]$V1
    bs$upr = b[, quantile(Length_µm, prob = c(0.975)), by = list(part,Morph)]$V1

    # mean +/-sd
    bs = b[, mean(Length_µm), by = list(part,Morph)]
    names(bs)[3] = 'median'
    bs$lwr = b[, mean(Length_µm)-sd(Length_µm), by = list(part,Morph)]$V1
    bs$upr = b[, mean(Length_µm)+sd(Length_µm), by = list(part,Morph)]$V1

    bs_f = bs[part == 'Flagellum']
    bs_f[,Midpiece_µm:= bs[part == 'Midpiece',.(median)]]

    bs_fm = bs[part == 'Midpiece']
    bs_fm[,Flagellum_µm:= bs[part == 'Flagellum',.(median)]]
    
    g0 =
    ggplot(h, aes(x = Flagellum_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
      geom_point(pch = 21)+
      geom_point(data = bs_f, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
      geom_segment(data = bs_f, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
      geom_segment(data = bs_fm, aes(x = Flagellum_µm, xend = Flagellum_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = bs_f, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_f, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      annotate(geom="text", x=80, y=29.5, label="Independent", size = 2, col = colors[1], hjust = 0) +
      annotate(geom="text", x=80, y=29, label="Satellite", size = 2, col = colors[2], hjust = 0) +
      annotate(geom="text", x=80, y=28.5, label="Faeder", size = 2, col = colors[3], hjust = 0) +
      annotate(geom="text", x=80, y=28, label="Mean +/-SD", size = 2, hjust = 0) +
      #scale_colour_viridis(discrete=TRUE)+
      ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

    bs_h = bs[part == 'Head']
    bs_h[,Midpiece_µm:= bs[part == 'Midpiece',.(median)]]

    bs_m = bs[part == 'Midpiece']
    bs_m[,Head_µm:= bs[part == 'Head',.(median)]]
    
    g1 =
    ggplot(h, aes(x = Head_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
      geom_point(pch = 21)+
      geom_point(data = bs_h, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
      geom_segment(data = bs_h, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
      geom_segment(data = bs_m, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = bs_h, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=25, y=29.5, label="Independent", size = 2, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=25, y=29, label="Satellite", size = 2, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=25, y=28.5, label="Faeder", size = 2, col = colors[3], hjust = 0) +
      #annotate(geom="text", x=25, y=28, label="Mean +/-SD", size = 2, hjust = 0) +
      #scale_colour_viridis(discrete=TRUE)+
      ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

    bs_h_ = bs[part == 'Head']
    bs_h_[,Flagellum_µm:= bs[part == 'Flagellum',.(median)]]

    bs_f = bs[part == 'Flagellum']
    bs_f[,Head_µm:= bs[part == 'Head',.(median)]]
      
    g2 =
    ggplot(h, aes(x = Head_µm, y = Flagellum_µm, col = Morph, fill = Morph)) +
      geom_point(pch = 21)+
      geom_point(pch = 21)+
      geom_point(data = bs_h_, aes(x = median, y = Flagellum_µm), col = 'black', size = 2) +
      geom_segment(data = bs_h_, aes(x = lwr, xend = upr, y =Flagellum_µm, yend =Flagellum_µm), col = "black" ) +
      geom_segment(data = bs_f, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = bs_h_, aes(x = median, y = Flagellum_µm, col = Morph), size = 1) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h_, aes(x = median, y = Flagellum_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #scale_colour_viridis(discrete=TRUE)+
      theme_bw() +
      theme(legend.position = "none",
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )   

    bs_n = bs[part == 'Nucleus']
    bs_n[,Acrosome_µm:= bs[part == 'Acrosome',.(median)]]

    bs_a = bs[part == 'Acrosome']
    bs_a[,Nucleus_µm:= bs[part == 'Nucleus',.(median)]] 
    
    g3 =
    ggplot(h, aes(x = Nucleus_µm, y = Acrosome_µm, col = Morph, fill = Morph)) +
     geom_point(pch = 21)+
      geom_point(data = bs_n, aes(x = median, y = Acrosome_µm), col = 'black', size = 2) +
      geom_segment(data = bs_n, aes(x = lwr, xend = upr, y =Acrosome_µm, yend =Acrosome_µm), col = "black" ) +
      geom_segment(data = bs_a, aes(x = Nucleus_µm, xend = Nucleus_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = bs_n, aes(x = median, y = Acrosome_µm, col = Morph), size = 1) +
      # geom_point(data = bs_n, aes(x = median, y = Acrosome_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #scale_colour_viridis(discrete=TRUE)+
      theme_bw() +
      theme(legend.position = "none",
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )  
  # male means
    as = a[, quantile(Length_avg, prob = c(0.5)), by = list(part,Morph)]
    names(as)[3] = 'median'
    as$lwr = a[, quantile(Length_avg, prob = c(0.025)), by = list(part,Morph)]$V1
    as$upr = a[, quantile(Length_avg, prob = c(0.975)), by = list(part,Morph)]$V1

    # mean +/-sd
    as = a[, mean(Length_avg), by = list(part,Morph)]
    names(as)[3] = 'median'
    as$lwr = a[, mean(Length_avg)-sd(Length_avg), by = list(part,Morph)]$V1
    as$upr = a[, mean(Length_avg)+sd(Length_avg), by = list(part,Morph)]$V1

    as_f = as[part == 'Flagellum']
    as_f[,Midpiece_µm:= as[part == 'Midpiece',.(median)]]

    as_fm = as[part == 'Midpiece']
    as_fm[,Flagellum_µm:= as[part == 'Flagellum',.(median)]]
    
    g0a =
    ggplot(ha, aes(x = Flagellum_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
      geom_point(pch = 21)+
      geom_point(data = as_f, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
      geom_segment(data = as_f, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
      geom_segment(data = as_fm, aes(x = Flagellum_µm, xend = Flagellum_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = as_f, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = aas_f, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
     #scale_colour_viridis(discrete=TRUE)+
      ggtitle('Male means') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

    as_h = as[part == 'Head']
    as_h[,Midpiece_µm:= as[part == 'Midpiece',.(median)]]

    as_m = as[part == 'Midpiece']
    as_m[,Head_µm:= as[part == 'Head',.(median)]]

    g1a =
    ggplot(ha, aes(x = Head_µm, y = Midpiece_µm, col = Morph, fill = Morph)) +
      geom_point(pch = 21)+
      geom_point(data = as_h, aes(x = median, y = Midpiece_µm), col = 'black', size = 2) +
      geom_segment(data = as_h, aes(x = lwr, xend = upr, y =Midpiece_µm, yend =Midpiece_µm), col = "black" ) +
      geom_segment(data = as_m, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = as_h, aes(x = median, y = Midpiece_µm, col = Morph), size = 1) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = as_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ggtitle('Male means') +
      theme_bw() +
      theme(legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        plot.title = element_text(size=9)
        ) 

    as_h_ = as[part == 'Head']
    as_h_[,Flagellum_µm:= as[part == 'Flagellum',.(median)]]

    as_f = as[part == 'Flagellum']
    as_f[,Head_µm:= as[part == 'Head',.(median)]]

    g2a =
    ggplot(ha, aes(x = Head_µm, y = Flagellum_µm, col = Morph, fill = Morph)) +
      geom_point(pch = 21)+
      geom_point(data = as_h_, aes(x = median, y = Flagellum_µm), col = 'black', size = 2) +
      geom_segment(data = as_h_, aes(x = lwr, xend = upr, y =Flagellum_µm, yend =Flagellum_µm), col = "black" ) +
      geom_segment(data = as_f, aes(x = Head_µm, xend = Head_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = as_h_, aes(x = median, y = Flagellum_µm, col = Morph), size = 1) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #scale_colour_viridis(discrete=TRUE)+
      theme_bw() +
      theme(legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        plot.title = element_text(size=9)
        )   

    as_n = as[part == 'Nucleus']
    as_n[,Acrosome_µm:= as[part == 'Acrosome',.(median)]]

    as_a = as[part == 'Acrosome']
    as_a[,Nucleus_µm:= as[part == 'Nucleus',.(median)]] 
    
    g3a =
    ggplot(ha, aes(x = Nucleus_µm, y = Acrosome_µm, col = Morph, fill = Morph)) +
      geom_point(pch = 21)+
      geom_point(data = as_n, aes(x = median, y = Acrosome_µm), col = 'black', size = 2) +
      geom_segment(data = as_n, aes(x = lwr, xend = upr, y =Acrosome_µm, yend =Acrosome_µm), col = "black" ) +
      geom_segment(data = as_a, aes(x = Nucleus_µm, xend = Nucleus_µm, y =lwr, yend =upr), col = "black" )+
      geom_point(data = as_n, aes(x = median, y = Acrosome_µm, col = Morph), size = 1) +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = as_n, aes(x = median, y = Acrosome_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      theme_bw() +
      theme(legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        plot.title = element_text(size=9)
        ) 
  # export
    gg0 <- ggplotGrob(g0)
    gg1 <- ggplotGrob(g1)
    gg2 <- ggplotGrob(g2) 
    gg3 <- ggplotGrob(g3) 
    gg0a <- ggplotGrob(g0a)
    gg1a <- ggplotGrob(g1a)
    gg2a <- ggplotGrob(g2a) 
    gg3a <- ggplotGrob(g3a) 

    grid.draw(cbind(rbind(gg0,gg1,gg2, gg3, size = "first"), rbind(gg0a,gg1a,gg2a, gg3a, size = "first"), size = "first"))

      #grid.arrange(g1,g2)
    
    ggsave(here::here('Output/morpho_corWithMeans.png'),cbind(rbind(gg0,gg1,gg2, gg3, size = "first"), rbind(gg0a,gg1a,gg2a, gg3a, size = "first"), size = "first"), width = 10, height =20, units = 'cm') 
    bw[,morph123 :=ifelse(Morph == 'Independent', 1, ifelse(Morph == 'Satellite', 2,3))]             

    p = cloud(Tail ~ Head * Midpiece, data=bw, group = Morph, col = colors,
              key = list(text = list(c('Independent', 'Satellite','Faeder'),col=colors,cex=0.6)),
              pch = 16, cex = 0.8, alpha = 0.75, 
              xlab=list(cex=0.7),
              ylab=list(cex=0.7),
              zlab = list(cex=0.7) 
              )

    npanel <- c(4, 2)
    rotx <- c(-50, -80)
    rotz <- seq(30, 300, length = npanel[1]+1)  
    update(p[rep(1, prod(npanel))], layout = npanel,
           panel = function(..., screen) {
              panel.cloud(..., screen = list(z = rotz[current.column()],
                                                 x = rotx[current.row()]))
          })

    png(here::here('Output/morpho_cor_3D.png'),width = 20, height = 10, units = "cm", res = 600)
      update(p[rep(1, prod(npanel))], layout = npanel,
       panel = function(..., screen) {
        panel.cloud(..., screen = list(z = rotz[current.column()],
                x = rotx[current.row()]))
            })
       dev.off()

Models

Effect sizes

    # morpho  
      # for averages
        l = list()
        lp =list()
        lpr = list()
        lcv = list()
        lpcv = list()
        
        for(i in unique(a$part)){
          #i ='Nucleus'
          m = lm(scale(Length_avg) ~ Morph, a[part == i])
          #summary(m)
          #plot(allEffects(m))
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
          l[[i]]=data.frame(response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

          # get predictions
          m = lm(Length_avg ~ Morph, a[part == i])
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          newD=data.frame(Morph = unique(a$Morph)) # values to predict for
          X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here 
          newD$Length_avg <-(X%*%v) 
          predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
          for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                          predmatrix[predmatrix < 0] <- 0
                          newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                          newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                          #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
          newD$part=i
          lp[[i]] = data.table(newD)

          print(i)     
          }          
        for(i in unique(ar$part)){
          if(i == 'Midpiece'){ii = 'Midpiece_rel'}
          if(i == 'Flagellum'){ii = 'Flagellum_rel'}
          m = lm(scale(Length_rel) ~ Morph, ar[part == i])
          #summary(m)
          #plot(allEffects(m))
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
          l[[ii]]=data.frame(response=ii,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

          # get predictions
          m = lm(Length_rel ~ Morph, ar[part == i])
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          newD=data.frame(Morph = unique(a$Morph)) # values to predict for
          X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here 
          newD$Length_rel <-(X%*%v) 
          predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
          for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                          predmatrix[predmatrix < 0] <- 0
                          newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                          newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                          #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
          newD$part=i
          lpr[[i]] = data.table(newD)

          print(i)     
          }          
        for(i in unique(cv_$part)){
          #i ='Nucleus'
          m = lm(scale(CV) ~ Morph, cv_[part == i])
          #summary(m)
          #plot(allEffects(m))
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
          lcv[[i]]=data.frame(response=paste('CV',i),effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

          # get predictions
          m = lm(CV ~ Morph, cv_[part == i])
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@coef, 2, quantile, prob=c(0.5))
          newD=data.frame(Morph = unique(a$Morph)) # values to predict for
          X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here 
          newD$Length_avg <-(X%*%v) 
          predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
          for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                          predmatrix[predmatrix < 0] <- 0
                          newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                          newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                          #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
          newD$part=i
          lpcv[[i]] = data.table(newD)

          print(i)     
          }          
        
        ll = data.table(do.call(rbind,l) ) 
        ll[, response := factor(response, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))] 
        ll[, effect := factor(effect, levels=rev(c("(Intercept)", "MorphSatellite", "MorphFaeder")))]
        ll[, unit := 'male average']

        llcv = data.table(do.call(rbind,lcv) ) 
        llcv[, response := factor(response, levels=rev(c("CV Acrosome", "CV Nucleus", "CV Head", "CV Midpiece","CV Tail","CV Flagellum","CV Total")))] 
        llcv[, effect := factor(effect, levels=rev(c("(Intercept)", "MorphSatellite", "MorphFaeder")))]
        
        llp = data.table(do.call(rbind,lp) ) 
        llp[, part := factor(part, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total")))] 
        llp[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]

        llpr = data.table(do.call(rbind,lpr) ) 
        llpr[, part := factor(part, levels=rev(c("Midpiece","Flagellum")))] 
        llpr[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]

      # for single values 
        ls = list()
        lps = list()
        lpsr = list()
        for(i in unique(a$part)){
          #i ='Nucleus'
          m = lmer(scale(Length_µm) ~ Morph + (1|bird_ID), b[part == i])
          #summary(m)
          #plot(allEffects(m))
          bsim = sim(m, n.sim=5000) 
          v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 
          ls[[i]]=data.frame(response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
         
          # get predictions
          m = lmer(Length_µm ~ Morph + (1|bird_ID), b[part == i])
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
          newD=data.frame(Morph = unique(a$Morph)) # values to predict for
          X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here 
          newD$Length_µm <-(X%*%v) 
          predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
          for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@fixef[j,])}
                          predmatrix[predmatrix < 0] <- 0
                          newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                          newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                          #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
          newD$part=i
          lps[[i]] = data.table(newD)

          print(i)         
          }  
        for(i in unique(ar$part)){
          if(i == 'Midpiece'){ii = 'Midpiece_rel'}
          if(i == 'Flagellum'){ii = 'Flagellum_rel'}
          m = lmer(scale(Length_rel) ~ Morph + (1|bird_ID), br[part == i])
          #summary(m)
          #plot(allEffects(m))
          bsim = sim(m, n.sim=5000) 
          v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
          ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 
          ls[[ii]]=data.frame(response=ii,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])
         
          # get predictions
          m = lmer(Length_rel ~ Morph + (1|bird_ID), br[part == i])
          bsim = sim(m, n.sim=nsim) 
          v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
          newD=data.frame(Morph = unique(ar$Morph)) # values to predict for
          X <- model.matrix(~ Morph,data=newD) # exactly the model which was used has to be specified here 
          newD$Length_rel <-(X%*%v) 
          predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
          for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@fixef[j,])}
                          predmatrix[predmatrix < 0] <- 0
                          newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                          newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                          #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
          newD$part=i
          lpsr[[i]] = data.table(newD)

          print(i)     
          }          
        
        lls = data.table(do.call(rbind,ls) ) 
        lls[, response := factor(response, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))] 
        lls[, effect := factor(effect, levels=rev(c("(Intercept)", "MorphSatellite", "MorphFaeder")))]
        lls[, unit := 'single sperm']

        llps = data.table(do.call(rbind,lps) ) 
        llps[, part := factor(part, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))] 
        llps[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]

        llpsr = data.table(do.call(rbind,lpsr) ) 
        llpsr[, part := factor(part, levels=rev(c("Midpiece","Flagellum")))] 
        llpsr[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]

    # motility
      #chart.Correlation(d[, c('VAP', 'VSL','VCL', 'motileCount','motileCount_ln', 'NumberFields')], histogram=TRUE, pch=19)
      #mtext("Single sperm", side=3, line=3)
      lv = list()
      lvp =list()
      d[, motileCount_ln:=scale(log(motileCount))]
      dd = d[!Morph%in%'Zebra finch']
      #dd = a[part =='Acrosome']
      #dd[, motileCount_ln:=scale(log(motileCount))]
      # VAP
        m = lm(scale(VAP) ~ scale(log(motileCount)) + Morph, dd)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
        lv[['VAP']]=data.frame(response='Average-path (VAP)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lm(VAP ~ motileCount_ln + Morph, dd)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(dd$motileCount_ln), Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Average-path (VAP)'
        lvp[['vap']] = data.table(newD)   
      # VSL
        m = lm(scale(VSL) ~ scale(log(motileCount))  + Morph, dd)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
        lv[['VSL']]=data.frame(response='Straight-line (VSL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lm(VSL ~ motileCount_ln + Morph, dd)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(dd$motileCount_ln),Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Straight-line (VSL)'
        lvp[['VSL']] = data.table(newD)
      # VCL
        m = lm(scale(VCL) ~scale(log(motileCount)) + Morph, dd)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
        lv[['VCL']]=data.frame(response='Curvilinear (VCL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lm(VCL ~ motileCount_ln + Morph, dd)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(dd$motileCount_ln), Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Curvilinear (VCL)'
        lvp[['VCL']] = data.table(newD) 
                 
      llv = data.table(do.call(rbind,lv) ) 
      llv = llv[effect != 'scale(log(motileCount))' ]
      llv[, response := factor(response, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 
      llv[effect == '(Intercept)', effect:='Independent\n(Intercept)']
      llv[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
      llv[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
      
      llv[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))] 

      llvp = data.table(do.call(rbind,lvp) ) 
      llvp[, motility := factor(motility, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 
      llvp[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]
    # motility 2 - only one measure per bird
      #chart.Correlation(d[, c('VAP', 'VSL','VCL', 'motileCount','motileCount_ln', 'NumberFields')], histogram=TRUE, pch=19)
      #mtext("Single sperm", side=3, line=3)
      lvx = list()
      lvpx =list()
      d[, motileCount_ln:=scale(log(motileCount))]
      dd = d[!Morph%in%'Zebra finch']
      dd1 = dd[month == 'June']
      dd2 = dd[month == 'May']
      ddx = rbind(dd1,dd2[!bird_ID%in%dd1$bird_ID])
      # VAP
        m = lm(scale(VAP) ~ scale(log(motileCount)) + Morph, ddx)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
        lvx[['VAP']]=data.frame(response='Average-path (VAP)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lm(VAP ~ motileCount_ln + Morph, ddx)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(ddx$motileCount_ln), Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Average-path (VAP)'
        lvpx[['vap']] = data.table(newD)   
      # VSL
        m = lm(scale(VSL) ~ scale(log(motileCount))  + Morph, ddx)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
        lvx[['VSL']]=data.frame(response='Straight-line (VSL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lm(VSL ~ motileCount_ln + Morph, ddx)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(ddx$motileCount_ln),Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Straight-line (VSL)'
        lvpx[['VSL']] = data.table(newD)
      # VCL
        m = lm(scale(VCL) ~scale(log(motileCount)) + Morph, ddx)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
        lvx[['VCL']]=data.frame(response='Curvilinear (VCL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lm(VCL ~ motileCount_ln + Morph, ddx)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@coef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(ddx$motileCount_ln), Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Curvilinear (VCL)'
        lvpx[['VCL']] = data.table(newD) 
                 
      llvx = data.table(do.call(rbind,lvx) ) 
      llvx = llvx[effect != 'scale(log(motileCount))' ]
      llvx[, response := factor(response, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 
      llvx[effect == '(Intercept)', effect:='Independent\n(Intercept)']
      llvx[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
      llvx[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
      llvx[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))] 

      llvpx = data.table(do.call(rbind,lvpx) ) 
      llvpx[, motility := factor(motility, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 
      llvpx[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]
    # motility 3 - mixed model
      #chart.Correlation(d[, c('VAP', 'VSL','VCL', 'motileCount','motileCount_ln', 'NumberFields')], histogram=TRUE, pch=19)
      #mtext("Single sperm", side=3, line=3)
      lvq = list()
      lvpq =list()
      d[, motileCount_ln:=scale(log(motileCount))]
      dd = d[!Morph%in%'Zebra finch']
      # VAP
        m = lmer(scale(VAP) ~ scale(log(motileCount)) + Morph + (1|bird_ID), dd)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 
        lvq[['VAP']]=data.frame(response='Average-path (VAP)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lmer(VAP ~ motileCount_ln + Morph+ (1|bird_ID), dd)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(dd$motileCount_ln), Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@fixef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Average-path (VAP)'
        lvpq[['vap']] = data.table(newD)   
      # VSL
        m = lmer(scale(VSL) ~ scale(log(motileCount))  + Morph + (1|bird_ID), dd)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 
        lvq[['VSL']]=data.frame(response='Straight-line (VSL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lmer(VSL ~ motileCount_ln + Morph + (1|bird_ID), dd)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(dd$motileCount_ln),Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@fixef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Straight-line (VSL)'
        lvpq[['VSL']] = data.table(newD)
      # VCL
        m = lmer(scale(VCL) ~scale(log(motileCount)) + Morph + (1|bird_ID), dd)
        #summary(m)
        #plot(allEffects(m))
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
        ci = apply(bsim@fixef, 2, quantile, prob=c(0.025,0.975)) 
        lvq[['VCL']]=data.frame(response='Curvilinear (VCL)',effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

        # get predictions
        m = lmer(VCL ~ motileCount_ln + Morph + (1|bird_ID), dd)
        bsim = sim(m, n.sim=nsim) 
        v = apply(bsim@fixef, 2, quantile, prob=c(0.5))
        newD=data.frame(motileCount_ln = mean(dd$motileCount_ln), Morph = unique(b$Morph)) # values to predict for
        X <- model.matrix(~ motileCount_ln + Morph,data=newD) # exactly the model which was used has to be specified here 
        newD$pred <-(X%*%v) 
        predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
        for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@fixef[j,])}
                    predmatrix[predmatrix < 0] <- 0
                    newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                    newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                    #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
        newD$motility = 'Curvilinear (VCL)'
        lvpq[['VCL']] = data.table(newD) 
                 
      llvq = data.table(do.call(rbind,lvq) ) 
      llvq = llvq[effect != 'scale(log(motileCount))' ]
      llvq[, response := factor(response, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 
      llvq[effect == '(Intercept)', effect:='Independent\n(Intercept)']
      llvq[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
      llvq[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
      llvq[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))] 

      llvpq = data.table(do.call(rbind,lvpq) ) 
      llvpq[, motility := factor(motility, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 
      llvpq[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]
        llll = rbind(ll,lls)  
        llll[, unit := factor(unit, levels=rev(c("single sperm", "male average")))] 
        llll[effect == '(Intercept)', effect:='Independent\n(Intercept)']
        llll[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
        llll[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
        llll[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))] 
        
        g = 
        ggplot(llll, aes(y = effect, x = estimate, col = response, shape = unit)) +
          geom_vline(xintercept = 0, col = "grey30", lty =3)+
          geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.4) ) +
          #ggtitle ("Sim based")+
          geom_point(position = position_dodge(width = 0.4)) +
          #scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          #scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE))  +
          scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) + 
          scale_shape(guide = guide_legend(reverse = TRUE)) + 
          #scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
          labs(y = NULL ,x = "Standardized effect size",title = 'Sperm-part specific models') +
          #ylim(c(0,100))+
          #coord_flip()+
          theme_bw() +
          theme( #legend.position ="right",
                plot.title = element_text(size=7),
                legend.title=element_text(size=7), 
                legend.text=element_text(size=6),
                ##legend.spacing.y = unit(0.1, 'cm'), 
                legend.key.height= unit(0.5,"line"),
                #plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
                panel.grid = element_blank(),
                panel.border = element_blank(),
                panel.background = element_blank(),
                axis.line = element_line(colour = ax_lines, size = 0.25),
                axis.line.y = element_blank(),
                axis.ticks.y = element_blank(),
                axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
                axis.ticks.length = unit(1, "pt"),
                axis.text.x = element_text(colour="black", size = 7),
                axis.text.y=element_text(colour="black", size = 7),
                axis.title=element_text(size=9)
                )
        g

        ggsave(here::here('Output/morpho_effectSizes_virid.png'),g, width = 12, height =10, units = 'cm')
        #ggsave('Output/morpho_effectSizes_pair.png',g, width = 12, height =10, units = 'cm')

!!! Satellites seem to have the longest tails (and flagellums), faeders the longest midpieces, which translates into the differences in relative measurements. !!!


        llcv[effect == '(Intercept)', effect:='Independent\n(Intercept)']
        llcv[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
        llcv[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
        llcv[, effect := factor(effect, levels=c("Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))] 
        
        g = 
        ggplot(llcv, aes(y = effect, x = estimate, col = response)) +
          geom_vline(xintercept = 0, col = "grey30", lty =3)+
          geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.4) ) +
          #ggtitle ("Sim based")+
          geom_point(position = position_dodge(width = 0.4)) +
          #scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          #scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE))  +
          scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) + 
          #scale_shape(guide = guide_legend(reverse = TRUE)) + 
          #scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
          labs(y = NULL ,x = "Standardized effect size",title = 'Sperm-part specific models on male CV') +
          #ylim(c(0,100))+
          #coord_flip()+
          theme_bw() +
          theme( #legend.position ="right",
                plot.title = element_text(size=7),
                legend.title=element_text(size=7), 
                legend.text=element_text(size=6),
                ##legend.spacing.y = unit(0.1, 'cm'), 
                legend.key.height= unit(0.5,"line"),
                #plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
                panel.grid = element_blank(),
                panel.border = element_blank(),
                panel.background = element_blank(),
                axis.line = element_line(colour = ax_lines, size = 0.25),
                axis.line.y = element_blank(),
                axis.ticks.y = element_blank(),
                axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
                axis.ticks.length = unit(1, "pt"),
                axis.text.x = element_text(colour="black", size = 7),
                axis.text.y=element_text(colour="black", size = 7),
                axis.title=element_text(size=9)
                )
        g

        ggsave(here::here('Output/morpho_effectSizes_virid_CV.png'),g, width = 12, height =10, units = 'cm')
        g2 = 
        ggplot(llvx, aes(y = effect, x = estimate, col = response)) +
          geom_vline(xintercept = 0, col = "grey30", lty =3)+
          geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.4) ) +
          #ggtitle ("Sim based")+
          geom_point(position = position_dodge(width = 0.4)) +
          #scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          #scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE))  +
          scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) + 
          #scale_shape(guide = guide_legend(reverse = TRUE)) + 
          #scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
          labs(y = NULL ,x = "Standardized effect size",title = 'Velocity-type specific models on male velocity\nusing single recording per male (for all but three males from June)') +
          #ylim(c(0,100))+
          #coord_flip()+
          theme_bw() +
          theme( #legend.position ="right",
                plot.title = element_text(size=7),
                legend.title=element_text(size=7), 
                legend.text=element_text(size=6),
                ##legend.spacing.y = unit(0.1, 'cm'), 
                legend.key.height= unit(0.5,"line"),
                #plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
                panel.grid = element_blank(),
                panel.border = element_blank(),
                panel.background = element_blank(),
                axis.line = element_line(colour = ax_lines, size = 0.25),
                axis.line.y = element_blank(),
                axis.ticks.y = element_blank(),
                axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
                axis.ticks.length = unit(1, "pt"),
                axis.text.x = element_text(colour="black", size = 7),
                axis.text.y=element_text(colour="black", size = 7),
                axis.title=element_text(size=9)
                )
        g2

        ggsave(here::here('Output/motility_effectSizes_virid.png'),g2, width = 10, height =8, units = 'cm')
        g2 = 
        ggplot(llvq, aes(y = effect, x = estimate, col = response)) +
          geom_vline(xintercept = 0, col = "grey30", lty =3)+
          geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.4) ) +
          #ggtitle ("Sim based")+
          geom_point(position = position_dodge(width = 0.4)) +
          #scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          #scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE))  +
          scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) + 
          #scale_shape(guide = guide_legend(reverse = TRUE)) + 
          #scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
          labs(y = NULL ,x = "Standardized effect size",title = 'Velocity-type specific models on male velocity\nusing all recording and bird_ID as random intercept') +
          #ylim(c(0,100))+
          #coord_flip()+
          theme_bw() +
          theme( #legend.position ="right",
                plot.title = element_text(size=7),
                legend.title=element_text(size=7), 
                legend.text=element_text(size=6),
                ##legend.spacing.y = unit(0.1, 'cm'), 
                legend.key.height= unit(0.5,"line"),
                #plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
                panel.grid = element_blank(),
                panel.border = element_blank(),
                panel.background = element_blank(),
                axis.line = element_line(colour = ax_lines, size = 0.25),
                axis.line.y = element_blank(),
                axis.ticks.y = element_blank(),
                axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
                axis.ticks.length = unit(1, "pt"),
                axis.text.x = element_text(colour="black", size = 7),
                axis.text.y=element_text(colour="black", size = 7),
                axis.title=element_text(size=9)
                )
        g2

        ggsave(here::here('Output/motility_effectSizes_virid_mixed.png'),g2, width = 10, height =8, units = 'cm')

!!! Against prediction about selection for the fastest sperm, faeders might have most variable & slowest sperm !!!


Predictions with raw data

      g1 = 
          ggplot(a, aes(x = Morph, y = Length_avg)) +
          geom_dotplot(binaxis = 'y', stackdir = 'center',
                        position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1.5)+
          geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) + 
          geom_errorbar(data = llp, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
          geom_point(data = llp, aes(x = Morph, y =Length_avg), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
          #stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
          scale_color_viridis(discrete=TRUE)+
          facet_wrap(~part, scales = 'free_y', nrow = 2)+
          labs(title = 'Predictions from sperm part specific models\non male means') +
          ylab('Length [µm]') +
          theme_bw() +
          guides(x =  guide_axis(angle = -45)) +
          theme(legend.position = "none",
            axis.title.x = element_blank(), 
            #axis.text.x = element_blank(),
            plot.title = element_text(size=9)
            ) #panel.spacing.y = unit(0, "mm")) #,    
 
      g2 = 
          ggplot(b, aes(x = Morph, y = Length_µm)) +
          geom_dotplot(binaxis = 'y', stackdir = 'center',
                        position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1)+
          geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) + 
          geom_errorbar(data = llps, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
          geom_point(data = llps, aes(x = Morph, y =Length_µm), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
          #stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
          scale_color_viridis(discrete=TRUE)+
          facet_wrap(~part, scales = 'free_y', nrow = 2)+
          labs(title = 'Predictions from sperm part specific models\non single sperm') +
          ylab('Length [µm]') +
          theme_bw() +
          guides(x =  guide_axis(angle = -45)) +
          theme(legend.position = "none",
            axis.title.x = element_blank(), 
            axis.text.x = element_blank(),
            plot.title = element_text(size=9)
            ) #panel.spacing.y = unit(0, "mm")) #,     

      grid.draw(rbind(ggplotGrob(g2), ggplotGrob(g1), size = "last"))

      #grid.arrange(g1,g2)
      gg1 <- ggplotGrob(g1)
      gg2 <- ggplotGrob(g2) 
      ggsave(here::here('Output/morpho_predictions+boxPlots.png'),rbind(gg2,gg1, size = "last"), width = 15, height =15, units = 'cm')    
      g1 = 
          ggplot(br, aes(x = Morph, y = Length_rel)) +
          geom_dotplot(binaxis = 'y', stackdir = 'center',
                        position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1)+
          geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) + 
          geom_errorbar(data = llpsr, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
          geom_point(data = llpsr, aes(x = Morph, y =Length_rel), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
          #stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
          scale_color_viridis(discrete=TRUE)+
          facet_wrap(~part, scales = 'free_y', nrow = 2)+
          ylab('Relative length to total sperm length') +
          labs(title = 'Predictions from sperm part specific models\non single sperm') +
          theme_bw() +
          guides(x =  guide_axis(angle = -45)) +
          theme(legend.position = "none",
            plot.title = element_text(size=9)
            ) #panel.spacing.y = unit(0, "mm")) #,     

      g2 = 
          ggplot(ar, aes(x = Morph, y = Length_rel)) +
          geom_dotplot(binaxis = 'y', stackdir = 'center',
                        position = position_dodge(), col = 'grey', fill = 'lightgrey', dotsize = 1.5)+
          geom_boxplot(aes(col = Morph), fill = NA, alpha = 0.2) + 
          geom_errorbar(data = llpr, aes(ymin = lwr, ymax = upr), width = 0, position = position_dodge(width = 0.25), col = 'red' ) +
          geom_point(data = llpr, aes(x = Morph, y =Length_rel), position = position_dodge(width = 0.25), col = 'red', size = 0.75) +
          #stat_summary(fun.y=mean, geom="point", color="red", fill="red") +
          scale_color_viridis(discrete=TRUE)+
          facet_wrap(~part, scales = 'free_y', nrow = 2)+
          labs(title = '\non male means') +
          theme_bw() +
          guides(x =  guide_axis(angle = -45)) +
          theme(legend.position = "none",
            axis.title.y = element_blank(), 
            axis.text.y = element_blank(),
            plot.title = element_text(size=9)
            ) #panel.spacing.y = unit(0, "mm")) #,    
 
      gg1 <- ggplotGrob(g1)
      gg2 <- ggplotGrob(g2) 

      grid.draw(cbind(gg1, gg2, size = "last"))

      #grid.arrange(g1,g2)
      
      ggsave(here::here('Output/morpho_predictions+boxPlots_rel.png'),rbind(gg1,gg2, size = "last"), width = 7.5, height =15/4, units = 'cm')   

Does morpho predict motility?

    # not controlled for motile count  
      l = list()
      lp =list()
      # VAP
         for(i in unique(a$part)){
            #i ='Nucleus'
            ai = a[part == i]
            m = lm(scale(VAP) ~ Morph+scale(Length_avg), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            l[[paste('VAP',i)]]=data.frame(mot = 'VAP', response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VAP ~ Morph+Length_avg, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(Morph = j, Length_avg = seq(min(ai$Length_avg), max(ai$Length_avg), length.out = 200)) 
              }
            newD=do.call(rbind,nd)

            # values to predict for
            X <- model.matrix(~ Morph+Length_avg,data=newD) # exactly the model which was used has to be specified here 
            newD$VAP <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=i
            newD$mot = 'VAP'
            setnames(newD, old = 'VAP', new = 'motility')
            lp[[paste(i,newD$mot[1])]] = data.table(newD)

            print(paste(i,newD$mot[1]))     
            }          
      # VSL
         for(i in unique(a$part)){
            #i ='Nucleus'
            ai = a[part == i]
            m = lm(scale(VSL) ~ Morph+scale(Length_avg), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            l[[paste('VSL',i)]]=data.frame(mot = 'VSL', response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VSL ~ Morph+Length_avg, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(Morph = j, Length_avg = seq(min(ai$Length_avg), max(ai$Length_avg), length.out = 200)) 
              }
            newD=do.call(rbind,nd)
            X <- model.matrix(~ Morph+Length_avg,data=newD) # exactly the model which was used has to be specified here 
            newD$VSL <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=i
            newD$mot = 'VSL'
            setnames(newD, old = 'VSL', new = 'motility')
            lp[[paste(i,newD$mot[1])]] = data.table(newD)

            print(paste(i,newD$mot[1]))     
            }          
      # VCL
         for(i in unique(a$part)){
            #i ='Nucleus'
            ai = a[part == i]
            m = lm(scale(VCL) ~ Morph+scale(Length_avg), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            l[[paste('VCL',i)]]=data.frame(mot = 'VCL', response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VCL ~ Morph+Length_avg, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(Morph = j, Length_avg = seq(min(ai$Length_avg), max(ai$Length_avg), length.out = 200)) 
              }
            newD=do.call(rbind,nd)
            X <- model.matrix(~ Morph+Length_avg,data=newD) # exactly the model which was used has to be specified here 
            newD$VCL <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=i
            newD$mot = 'VCL'
            setnames(newD, old = 'VCL', new = 'motility')
            lp[[paste(i,newD$mot[1])]] = data.table(newD)

            print(paste(i,newD$mot[1]))     
            }          
                 
      ll = data.table(do.call(rbind,l) ) 
      ll[, response := factor(response, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))] 
      ll[effect == '(Intercept)', effect:='Independent\n(Intercept)']
      ll[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
      ll[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
      ll[effect == 'scale(Length_avg)', effect := 'Length_µm']
      
      ll[, effect := factor(effect, levels=c('Length_µm',"Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))] 

      ll[mot=='VCL', mot:='Curvilinear (VCL)']
      ll[mot=='VSL', mot:='Straight-line (VSL)']
      ll[mot=='VAP', mot:='Average-path (VAP)']
      ll[, mot := factor(mot, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 

      llp = data.table(do.call(rbind,lp) ) 
      llp[, part := factor(part, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))] 
      llp[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))]  
    # controlled for motile count  
      lz = list()
      lpz =list()
      lprz =list()
      a[,motileCount_ln_z := scale(log(motileCount))]
      ar[,motileCount_ln_z := scale(log(motileCount))]
      # VAP
         for(i in unique(a$part)){
            #i ='Nucleus'
            ai = a[part == i]
            m = lm(scale(VAP) ~ motileCount_ln_z+ Morph+scale(Length_avg), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            lz[[paste('VAP',i)]]=data.frame(mot = 'VAP', response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VAP ~ motileCount_ln_z + Morph+Length_avg, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(motileCount_ln_z = mean(ai$motileCount_ln_z),Morph = j, Length_avg = seq(min(ai$Length_avg), max(ai$Length_avg), length.out = 200)) 
              }
            newD=do.call(rbind,nd)

            # values to predict for
            X <- model.matrix(~ motileCount_ln_z + Morph+Length_avg,data=newD) # exactly the model which was used has to be specified here 
            newD$VAP <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=i
            newD$mot = 'VAP'
            setnames(newD, old = 'VAP', new = 'motility')
            lpz[[paste(i,newD$mot[1])]] = data.table(newD)

            print(paste(i,newD$mot[1]))     
            }          
         for(i in unique(ar$part)){
            #i ='Nucleus'
            if(i == 'Midpiece'){ii = 'Midpiece_rel'}
            if(i == 'Flagellum'){ii = 'Flagellum_rel'}
            ai = ar[part == i]
            m = lm(scale(VAP) ~ motileCount_ln_z+ Morph+scale(Length_rel), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            lz[[paste('VAP',ii)]]=data.frame(mot = 'VAP', response=ii,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VAP ~ motileCount_ln_z + Morph+Length_rel, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(motileCount_ln_z = mean(ai$motileCount_ln_z),Morph = j, Length_rel = seq(min(ai$Length_rel), max(ai$Length_rel), length.out = 200)) 
              }
            newD=do.call(rbind,nd)

            # values to predict for
            X <- model.matrix(~ motileCount_ln_z + Morph+Length_rel,data=newD) # exactly the model which was used has to be specified here 
            newD$VAP <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=ii
            newD$mot = 'VAP'
            setnames(newD, old = 'VAP', new = 'motility')
            lpz[[paste(ii,newD$mot[1])]] = data.table(newD)

            print(paste(ii,newD$mot[1]))     
            }            
      # VSL
         for(i in unique(a$part)){
            #i ='Nucleus'
            ai = a[part == i]
            m = lm(scale(VSL) ~ motileCount_ln_z + Morph+scale(Length_avg), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            lz[[paste('VSL',i)]]=data.frame(mot = 'VSL', response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VSL ~ motileCount_ln_z + Morph+Length_avg, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(motileCount_ln_z = mean(ai$motileCount_ln_z), Morph = j, Length_avg = seq(min(ai$Length_avg), max(ai$Length_avg), length.out = 200)) 
              }
            newD=do.call(rbind,nd)
            X <- model.matrix(~ motileCount_ln_z + Morph+Length_avg,data=newD) # exactly the model which was used has to be specified here 
            newD$VSL <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=i
            newD$mot = 'VSL'
            setnames(newD, old = 'VSL', new = 'motility')
            lpz[[paste(i,newD$mot[1])]] = data.table(newD)

            print(paste(i,newD$mot[1]))     
            }          
         for(i in unique(ar$part)){
            #i ='Nucleus'
            if(i == 'Midpiece'){ii = 'Midpiece_rel'}
            if(i == 'Flagellum'){ii = 'Flagellum_rel'}
            ai = ar[part == i]
            m = lm(scale(VSL) ~ motileCount_ln_z+ Morph+scale(Length_rel), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            lz[[paste('VSL',ii)]]=data.frame(mot = 'VSL', response=ii,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VSL ~ motileCount_ln_z + Morph+Length_rel, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(motileCount_ln_z = mean(ai$motileCount_ln_z),Morph = j, Length_rel = seq(min(ai$Length_rel), max(ai$Length_rel), length.out = 200)) 
              }
            newD=do.call(rbind,nd)

            # values to predict for
            X <- model.matrix(~ motileCount_ln_z + Morph+Length_rel,data=newD) # exactly the model which was used has to be specified here 
            newD$VSL <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=ii
            newD$mot = 'VSL'
            setnames(newD, old = 'VSL', new = 'motility')
            lpz[[paste(ii,newD$mot[1])]] = data.table(newD)

            print(paste(ii,newD$mot[1]))     
            }            
      # VCL
         for(i in unique(a$part)){
            #i ='Nucleus'
            ai = a[part == i]
            m = lm(scale(VCL) ~ motileCount_ln_z + Morph+scale(Length_avg), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            lz[[paste('VCL',i)]]=data.frame(mot = 'VCL', response=i,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VCL ~ motileCount_ln_z + Morph+Length_avg, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(motileCount_ln_z = mean(ai$motileCount_ln_z), Morph = j, Length_avg = seq(min(ai$Length_avg), max(ai$Length_avg), length.out = 200)) 
              }
            newD=do.call(rbind,nd)
            X <- model.matrix(~ motileCount_ln_z + Morph+Length_avg,data=newD) # exactly the model which was used has to be specified here 
            newD$VCL <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=i
            newD$mot = 'VCL'
            setnames(newD, old = 'VCL', new = 'motility')
            lpz[[paste(i,newD$mot[1])]] = data.table(newD)

            print(paste(i,newD$mot[1]))     
            }          
         for(i in unique(ar$part)){
            #i ='Nucleus'
            if(i == 'Midpiece'){ii = 'Midpiece_rel'}
            if(i == 'Flagellum'){ii = 'Flagellum_rel'}
            ai = ar[part == i]
            m = lm(scale(VCL) ~ motileCount_ln_z+ Morph+scale(Length_rel), ai)
            #summary(m)
            #plot(allEffects(m))
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            ci = apply(bsim@coef, 2, quantile, prob=c(0.025,0.975)) 
            lz[[paste('VCL',ii)]]=data.frame(mot = 'VCL', response=ii,effect=rownames(coef(summary(m))),estimate=v, lwr=ci[1,], upr=ci[2,])

            # get predictions
            m = lm(VCL ~ motileCount_ln_z + Morph+Length_rel, ai)
            bsim = sim(m, n.sim=nsim) 
            v = apply(bsim@coef, 2, quantile, prob=c(0.5))
            nd = list()
            for(j in unique(ai$Morph)){
              #j=1
              nd[[j]] = data.frame(motileCount_ln_z = mean(ai$motileCount_ln_z),Morph = j, Length_rel = seq(min(ai$Length_rel), max(ai$Length_rel), length.out = 200)) 
              }
            newD=do.call(rbind,nd)

            # values to predict for
            X <- model.matrix(~ motileCount_ln_z + Morph+Length_rel,data=newD) # exactly the model which was used has to be specified here 
            newD$VCL <-(X%*%v) 
            predmatrix <- matrix(nrow=nrow(newD), ncol=nsim)
            for(j in 1:nsim) {predmatrix[,j] <- (X%*%bsim@coef[j,])}
                            predmatrix[predmatrix < 0] <- 0
                            newD$lwr <- apply(predmatrix, 1, quantile, prob=0.025)
                            newD$upr <- apply(predmatrix, 1, quantile, prob=0.975)
                            #newD$pred <- apply(predmatrix, 1, quantile, prob=0.5)
            newD$part=ii
            newD$mot = 'VCL'
            setnames(newD, old = 'VCL', new = 'motility')
            lpz[[paste(ii,newD$mot[1])]] = data.table(newD)

            print(paste(ii,newD$mot[1]))     
            }          
               
      llz = data.table(do.call(rbind,lz) ) 
      llz[, response := factor(response, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))] 
      llz[effect == '(Intercept)', effect:='Independent\n(Intercept)']
      llz[effect == 'MorphSatellite', effect := 'Satellite\n(relative to Independent)']
      llz[effect == 'MorphFaeder', effect := 'Faeder\n(relative to Independent)']
      llz[effect == 'scale(Length_avg)', effect := 'Length_µm']
      llz[effect == 'scale(Length_rel)', effect := 'Length_µm'] # dummy variable
      
      llz[, effect := factor(effect, levels=c('Length_µm','motileCount_ln_z',"Faeder\n(relative to Independent)","Satellite\n(relative to Independent)","Independent\n(Intercept)"))] 

      llz[mot=='VCL', mot:='Curvilinear (VCL)']
      llz[mot=='VSL', mot:='Straight-line (VSL)']
      llz[mot=='VAP', mot:='Average-path (VAP)']
      llz[, mot := factor(mot, levels=rev(c("Curvilinear (VCL)", "Straight-line (VSL)", "Average-path (VAP)")))] 

      llpz = data.table(do.call(rbind,lp) ) 
      llpz[, part := factor(part, levels=rev(c("Acrosome", "Nucleus", "Head", "Midpiece","Tail","Flagellum","Total","Midpiece_rel","Flagellum_rel")))] 
      llpz[, Morph := factor(Morph, levels=rev(c("Independent", "Satellite", "Faeder")))] 
        g10 = 
        ggplot(llz[effect%in%'Length_µm'], aes(y = response, x = estimate, col = response, shape = mot)) +
          geom_errorbar(aes(xmin = lwr, xmax = upr, col = response), width = 0.1, position = position_dodge(width = 0.5) ) +
          #ggtitle ("Sim based")+
          geom_point(position = position_dodge(width = 0.5)) +
          geom_vline(xintercept = 0, col = "red", lty =1)+
          #scale_colour_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          #scale_fill_brewer(type = 'qual', palette = 'Paired',guide = guide_legend(reverse = TRUE))+
          scale_color_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE))  +
          scale_fill_viridis(discrete=TRUE,guide = guide_legend(reverse = TRUE)) + 
          scale_shape(guide = guide_legend(reverse = TRUE, name = 'Motility')) + 
          #scale_x_continuous(limits = c(-2, 2), expand = c(0, 0), breaks = seq(-2,2, by = 1), labels = seq(-2,2, by = 1)) +
          labs(y = 'Predictor' , x = "Standardized effect size") +
          guides(color=FALSE, fill = FALSE) +
          #coord_fixed(ratio = 1.5)+
          #ylim(c(0,100))+
          #coord_flip()+
          theme_bw() +
          theme( #legend.position ="right",
                plot.title = element_text(size=7),
                legend.title=element_text(size=7), 
                legend.text=element_text(size=6),
                ##legend.spacing.y = unit(0.1, 'cm'), 
                legend.key.height= unit(0.5,"line"),
                #plot.margin = margin(b = 0.5, l = 0.5, t = 0.5, r =0.5, unit =  "pt"),
                panel.grid = element_blank(),
                panel.border = element_blank(),
                panel.background = element_blank(),
                axis.line = element_line(colour = ax_lines, size = 0.25),
                axis.line.y = element_blank(),
                axis.ticks.y = element_blank(),
                axis.ticks.x= element_line( colour = ax_lines, size = 0.25),
                axis.ticks.length = unit(1, "pt"),
                axis.text.x = element_text(colour="black", size = 7),
                axis.text.y=element_text(colour="black", size = 7),
                axis.title=element_text(size=9)
                )
          g10

          ggsave(here::here('Output/morpho_motil_effectSizes_virid_simple_controlled_MotileCount.png'),g10, width = 10, height =7, units = 'cm')
  ga1 = 
  ggplot(aw, aes(x = Acrosome, y = VAP)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #scale_colour_manual(values = colors)+
      #scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=3, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=3, y=13, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=3, y=11, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      ggtitle('Male means') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  
  ga2 = 
  ggplot(aw, aes(x = Acrosome, y = VSL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )    
  ga3 = 
  ggplot(aw, aes(x = Acrosome, y = VCL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
     
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

  gn1 = 
  ggplot(aw, aes(x = Nucleus, y = VAP)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
             #scale_colour_viridis(discrete=TRUE)+
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  gn2 = 
  ggplot(aw, aes(x = Nucleus, y = VSL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
            #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )    
  gn3 = 
  ggplot(aw, aes(x = Nucleus, y = VCL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

  gm1 = 
  ggplot(aw, aes(x = Midpiece, y = VAP)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
            #scale_colour_viridis(discrete=TRUE)+
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  gm2 = 
  ggplot(aw, aes(x = Midpiece, y = VSL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )    
  gm3 = 
  ggplot(aw, aes(x = Midpiece, y = VCL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  
  gt1 = 
  ggplot(aw, aes(x = Tail, y = VAP)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
            #scale_colour_viridis(discrete=TRUE)+
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  gt2 = 
  ggplot(aw, aes(x = Tail, y = VSL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
        #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )    
  gt3 = 
  ggplot(aw, aes(x = Tail, y = VCL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
       #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

  gh1 = 
  ggplot(aw, aes(x = Head, y = VAP)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  gh2 = 
  ggplot(aw, aes(x = Head, y = VSL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        plot.title = element_text(size=9)
        )    
  gh3 = 
  ggplot(aw, aes(x = Head, y = VCL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

  gf1 = 
  ggplot(aw, aes(x = Flagellum, y = VAP)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  gf2 = 
  ggplot(aw, aes(x = Flagellum, y = VSL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )    
  gf3 = 
  ggplot(aw, aes(x = Flagellum, y = VCL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

  go1 = 
  ggplot(aw, aes(x = Total, y = VAP)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
            #scale_colour_viridis(discrete=TRUE)+
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 
  go2 = 
  ggplot(aw, aes(x = Total, y = VSL)) +
      geom_smooth(method = 'rlm') +
      geom_point(pch = 21)+
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        )    
  go3 = 
  ggplot(aw, aes(x = Total, y = VCL)) +
      geom_point(pch = 21)+
      geom_smooth(method = 'rlm') +
      scale_colour_manual(values = colors)+
      scale_fill_manual(values = alpha(colors, 0.4))+
      #geom_point(data = bs_h, aes(x = median, y = Midpiece_µm),position = position_dodge(width = 0.25), col = 'black', size = 2) +
      #ylim(c(20,28)) +
      #xlim(c(24,36)) +
      #annotate(geom="text", x=28, y=15, label="Independent", size = 3, col = colors[1], hjust = 0) +
      #annotate(geom="text", x=28, y=14, label="Satellite", size = 3, col = colors[2], hjust = 0) +
      #annotate(geom="text", x=28, y=13, label="Faeder", size = 3, col = colors[3], hjust = 0) +
            #scale_colour_viridis(discrete=TRUE)+
      #ggtitle('Single measurements') +
      theme_bw() +
      theme(
        #legend.position=c(.9,.9),  
        legend.position = "none",
        axis.title.y = element_blank(), 
        axis.text.y = element_blank(),
        #axis.title.x = element_blank(), 
        #axis.text.x = element_blank(),
        plot.title = element_text(size=9)
        ) 

  grid.draw(cbind(
    rbind(ggplotGrob(ga1), ggplotGrob(ga2), ggplotGrob(ga3), size = "first"),
    rbind(ggplotGrob(gn1), ggplotGrob(gn2), ggplotGrob(gn3), size = "first"),
    rbind(ggplotGrob(gm1), ggplotGrob(gm2), ggplotGrob(gm3), size = "first"),
    rbind(ggplotGrob(gt1), ggplotGrob(gt2), ggplotGrob(gt3), size = "first"),
    rbind(ggplotGrob(gh1), ggplotGrob(gh2), ggplotGrob(gh3), size = "first"),
    rbind(ggplotGrob(gf1), ggplotGrob(gf2), ggplotGrob(gf3), size = "first"),
    rbind(ggplotGrob(go1), ggplotGrob(go2), ggplotGrob(go3), size = "first"),
    size = "first")
    )

  ggsave(here::here('Output/morpho_motil_simple.png'),cbind(
          rbind(ggplotGrob(ga1), ggplotGrob(ga2), ggplotGrob(ga3), size = "first"),
          rbind(ggplotGrob(gn1), ggplotGrob(gn2), ggplotGrob(gn3), size = "first"),
          rbind(ggplotGrob(gm1), ggplotGrob(gm2), ggplotGrob(gm3), size = "first"),
          rbind(ggplotGrob(gt1), ggplotGrob(gt2), ggplotGrob(gt3), size = "first"),
          rbind(ggplotGrob(gh1), ggplotGrob(gh2), ggplotGrob(gh3), size = "first"),
          rbind(ggplotGrob(gf1), ggplotGrob(gf2), ggplotGrob(gf3), size = "first"),
          rbind(ggplotGrob(go1), ggplotGrob(go2), ggplotGrob(go3), size = "first"),
          size = "first"), 
          width = 18, height =8, units = 'cm') 

!!! Perhaps a tendency for longer acrosomes making slower sperm and longer other traits making faster sperm !!!

# End